arXiv:l 508.05406V 1 [stat.ME] 21 Aug 2015 


A Note on Spherical Needlets 

Minjie Fan * 

Department of Statistics, University of California, Davis 


August 25, 2015 


Abstract 

Compared with the traditional spherical harmonics, the spherical needlets are a new genera¬ 
tion of spherical wavelets that possess several attractive properties. Their double localization in 
both spatial and frequency domains empowers them to easily and sparsely represent functions 
with small spatial scale features. This paper is divided into two parts. First, it reviews the 
spherical harmonics and discusses their limitations in representing functions with small spatial 
scale features. To overcome the limitations, it introduces the spherical needlets and their attrac¬ 
tive properties. In the second part of the paper, a Matlab package for the spherical needlets is 
presented. The properties of the spherical needlets are demonstrated by several examples using 
the package. 


1 Spherical Harmonics 

Plainly speaking, the spherical harmonics {T/m? / = 0,1, • • •, m = —•••,/} are a natural extension 
of the Fourier basis functions to the domain of the unit sphere For convenience, we employ 
(0, 0), 0 < 0 < TT, 0 < 0 < 27r, where 9 is the co-latitude and (j) is the longitude, as the spherical 
coordinate of x G The spherical harmonics share two crucial characteristics as the Fourier basis 
functions. First, they are eigenfunctions of the eigenvalue problem 


As2y = -Ay, 

where Ag 2 is the Laplace-Beltrami operator defined on i.e., 

A ^ 9 f . ^d\ 1 

“ sin e de V ™ ^de)^ sin^ 6 d(p‘^ ' 

Besides, they constitute an orthonormal basis of the Hilbert space Thus, for any function 

T in it can be uniquely expanded as 

oo I 

nx) = E E (3) 
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^im — / T{x^Yifyi{x)dx^ 

which are called spherical harmonic coefficients. 

The spherical harmonics can be represented in terms of the associated Legendre polynomials. 

Definition 1. For any / = 0, 1, • • •, m = —•••,/; the spherical harmonics 

Y,,M <P) = I ™ ^ “ (5) 

l(-l)"'r ro<0, 

where Pim represents the associated Legendre polynomial with subscripts I and m. 

The subscripts I and m entail the pattern of the waves along the latitude and longitude, where 
m is the longitudinal wave number and (/ — m + l)/2 is the latitudinal wave number. Proposition 
1 gives the addition formula for the spherical harmonics. 

Proposition 1. (Addition Formula) For any G 

^ _ 9/ -k 1 

Yim{x)Yim{y) = ( 6 ) 

m=—l 

where (•, •) denotes the inner product on and Pi represents the l-th Legendre polynomial. 


2 Spherical Needlets 

The spherical harmonics are globally supported on the sphere and thus suffer from the same dif¬ 
ficulties as the Fourier basis functions. For example, the Gibbs phenomenon may occur when a 
spherical function with small spatial scale features is approximated by a finite series of spherical 
harmonics. In this case, noticeable spurious global oscillations will be introduced and higher order 
of spherical harmonics are required to achieve a better approximation. 

In this paper, we shall review a new generation of spherical wavelets, called spherical needlets 
(Narcowich et al.^ 2006; Marinucci et al.^ 2008; Baldi et al.^ 2009; Marinucci and Peccati, 2011). 
They are not only exactly localized at a finite number of frequencies, but also decay quasi- 
exponentially fast away from their global maximum. 


2.1 Construction of Spherical Needlets 


The construction of the spherical needlets is based on two main ideas, which are the discretization 
of the sphere by an exact quadrature formula and a Littlewood-Paley decomposition. Theorem 1 
gives the exact quadrature formula. 


Theorem 1. Denote FLi as the space spanned by {T/m : m = —•••,/}; and let fCi — 

For any I G N, there exist a finite subset Xi — : A: = 1, • • •, n/} and positive real numbers 

{Xlk : A: = 1, • • •, ni} such that 


P ni 

/ f{x)dx 


( 7 ) 


for any / G /C;. Here $,ik and Xik are called cubature points and cubature weights, respectively. 
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The quadrature formula discretizes the sphere into cubature points and cubature weights, based 
on which, we present the definition of the spherical needlets. 

Definition 2. For a given frequency j G No and the corresponding cubature points ^jk and cubature 
weights Xjk (which are obtained by Theorem 1 with I = the spherical needlets with 

frequency j are defined as 


^jk ^ ^ ^ 


I 




m=—l 


( 8 ) 


= a/A^ T, (by Proposition 1), 

where x G S > 1 is a parameter and 6(-) is a window function satisfying 

1. b{-) > 0 in and it equals zero otherwise. 

2. For any ri>l, 

oo 

j=o 


(9) 


3. b{-) is M times continuously differentiable for some M = 1, 2, • • • or M = oo. 

Equation (8) implies that the spherical needlets are real-valued functions. The window function 
bf) plays the role of the Littlewood-Paley decomposition, i.e., it decomposes the frequency domain 
into several overlapping intervals {B^~^^ j = 0,1, • • •. Here represents the location (i.e., 

center) of while j determines to what extent f)jk is spatially localized. The larger j, the 
finer Varying the values of and j has the same effect as the translation and dilation in a 
multiresolution analysis. Thus, compared with the spherical harmonics, the spherical needlets are 
more suited to represent spherical functions with sharp local peaks or valleys. 


2.2 Properties of Spherical Needlets 

Compared with other spherical wavelets, the spherical needlets possess several attractive properties: 


1 . 


They are exactly localized in the frequency domain since the window function bf) has compact 
support. 


2. They are also localized in the spatial domain. Specifically, we have 


\'>Pjk{x)\ < 


_ cmB^ _ 

[1 + Bt arccos((^jfc, x))]^ 


for any ® € S^, 


( 10 ) 


where M G N such that 6 is M times continuously differentiable, and cm is a constant 
regardless of j and k. 

3. The spherical needlets (together with the first spherical harmonic Iqo = \/l/(47r)) constitute 
a Parseval tight frame, i.e., for any T G L^(S^), 


^liriii < + a^o < b\\t\\1 (11) 

j,k 
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where A = B = and 


{T,^jk) = [ T{x)ipjk{x)dx, 
which are called needlet coefficients and denoted as 


/ Tj]^(x)dx 



m (^jk) 



m 


( 12 ) 


where aim stre the spherical harmonic coefficients defined by equation (4). 

Equation (11) entails that (1) the needlet expansion preserves the “energy” of the function 
in the sense of L2-norm; (2) the dual frame coincides with the frame itself, and thus the 
spherical needlets have the perfect reconstruction property of orthonormal bases 


T{x) = aoo>oo(£c) + y^l3jk'>Pjk{x). 

j^k 


(13) 


4. The spherical needlets are almost orthogonal, i.e., if \j — j'\ > 2, 




'k'j 


/ 

JS2 


ipjk{x)'ipjik'{x)dx = 0 . 


5. The spherical needlets 'ijjjk and are asymptotically uncorrelated as the frequency increases 
and the distance between them remains fixed, i.e.. 


llV’jfcllllV’ifc'll 


<_ Cm _ 

“ [1 + arccos((^jfc, ijk'))]^ ’ 


where 

{i>jk,i>jk') = y ^jk^jk' T(si) 4^ Pi{{^jk,^jk'))- 

The proof is omitted here, and we refer to Lemma 3 in Baldi et al (2009). 


2.3 Implementation 

The needlet coefficients can be obtained by equation (12), which involves the spherical harmonic 
coefficients. When the sampling locations are on an equiangular grid, we can utilize the package 
S2Kit (Healy Jr et a/., 2003; Kostelec and Rockmore, 2004) to perform a fast spherical harmonic 
transform. When they are irregularly spaced, the spherical harmonic coefficients can be estimated 
empirically by 

N 

1=1 
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where i = 1, • • •, are the sampling locations, and wi is an appropriate cubature weight deter¬ 
mined by the surface area associated with xi. For example, we can estimate wi to be the area of 
the corresponding spherical polygon in the Voronoi diagram of all the sampling locations. 

It is assumed that the cubature points ^jk and cubature weights \jk are provided by the 
HEALPix discretization of the sphere (Gorski et al.^ 2005). Plainly speaking, the HEALPix grid 
discretizes the sphere into A^pix pixels with equal area, where A^pix = ^‘^^side ^side is required 
to be a power of two that measures the resolution of the discretization. We specify the cubature 
points as the center of the pixels, and the cubature weights as 

^jk — WT 7 (15) 


which is the equal area of the pixels. 

Suppose the highest frequency of the spherical harmonics that can be reliably extracted from 
the observations is /max- The maximum index jmax of Iho spherical needlets is the maximal j such 
that < /max- For each j between 0 and jmax (inclusively), the value of N^ide is determined by 

the inequality < 2A/side such that the quadrature formula in Theorem 1 holds approximately 

(Pietrobon et a/., 2010). There are other ways of discretizating the sphere with an exact quadrature 
formula and (almost) equal cubature weights, such as GLESP (Doroshkevich et a/., 2005) and 
(symmetric) spherical t-designs (Womersley, 2015). 

Computing the needlet coefficients by equation (12) is equivalent to performing an inverse 
spherical harmonic transform. The latter can be done in an efficient way by utilizing the distinct 
structure of the HEALPix grid. Specifically, the cubature points are arranged on a number of iso¬ 
latitude rings. The points on each ring share the same value of the co-latitude. Thus, the associated 
Legendre polynomials in the definition of the spherical harmonics only need to be evaluated once 
for all the points on each ring. Moreover, the grid is symmetric with respect to the equator, 
which reduces by half the computational cost of evaluating the associated Legendre polynomials. 
The inverse fast Fourier transform is also incorporated to further reduce the computational cost 
(see Appendix). Other computational techniques, such as parallelization and butterfly matrix 
compression are discussed in Hupca et al (2012); Seljebotn (2012). 

3 NeedMat: A Mat lab Package for Spherical Needlets 

Currently, there are several packages available for the spherical needlets, such as NeedATool 
(Pietrobon et a/., 2010) and S2LET (Leistedt et al.^ 2013). The former is written in Fortran 
and the latter is not easy to install because of its dependencies. In this section, we shall introduce 
a Matlab package for the spherical needlets, called NeedMat, which is easy to install and use. In 
the package, the window function 6(-) is constructed according to Marinucci et al (2008). 

3.1 Plot of Spherical Needlets 

The spherical needlet ^jk can be plotted by the function 

plot.needlets(B, j, k, res); 

where the argument res is a parameter controlling the resolution of the plots. Figures 1-2 show 
the plots of the spherical needlet with j = 3 and k = 100. 
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Figure 1: Plot of the spherical needlet 'ijjjk with j = 3 and k = 100 as a function of the great-circle 
distance between and x. Note that this plot remains the same for an arbitrary value of k. 
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Figure 2: Plot of the spherical needlet with j = 3 and k = 100 on the sphere. The sphere has 
been projected to an ellipse by the Hammer projection. 
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Figure 3: Computation time of the spherical needlet transform for different /max- 
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Figure 4: Observations of </3(-) on a perturbed HEALPix grid with 768 grid points. 
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3.2 Spherical Needlet Transform 

The spherical needlet transform can be performed by the function 
beta = spneedlet.tran(aim, l_max, B); 

We examine the computation time of the spherical needlet transform for different /max with aim 
randomly generated, as shown in Figure 3. All tests were run on a regular laptop with a 2.4 GHz 
Intel Core i5 processor. 

3.3 Example: Approximation of a Wendland Radial Basis Function 

In this subsection, we present an example to demonstrate the usage of NeedMat and the superiority 
of the spherical needlets to the spherical harmonics in capturing small spatial scale features. The 
function to be approximated is 


(p{x) = (po (arccos((|o,a^))//>), 

where is the center of the function, p is a spatial scale parameter, and is the Wendland 

radial basis function (Wendland, 1995) in 

1 — d)^{3b(f + 18d + 3)/3 for 0 < d < 1 



otherwise. 


We specify = (^/2,vr) in spherical coordinates and p — 7r/4. This specification implies that (p(-) 
has compact support, which is the spherical cap centered at with radius 7r/4. The observations 


are on a perturbed HEALPix grid with 768 grid points, as shown in Figure 4. The spherical 
harmonic coefficients are estimated by equation (14), where Wi are determined by the Voronoi 
diagram. It is achieved by performing 

aim = spharmonic_tran_irr(theta, phi, f_wend, l_max); 

where /max is set as 16. 

We reconstruct the function by the estimated spherical harmonic coefficients, as shown in Figure 
5. The reconstructed function has noticeable spurious global oscillations, which are the artifacts 
induced by the spherical harmonics. Based on the estimated spherical harmonic coefficients, we 
obtain the needlet coefficients by equation (12). It is achieved by performing 

beta = spneedlet_tran(aim, l_max, B); 

Figure 7 shows the estimated probability density functions of the needlet coefficients at frequencies 
j — 2,3,4. The spiky probability density functions entail that the spherical needlets can parsimo¬ 
niously represent this spatially localized function, which is one of the advantages of the spherical 
needlets over the spherical harmonics. Theoretically, without applying any shrinking operator to 
the needlet coefficients, the function reconstructed by the spherical needlets should be exactly the 
same as that reconstructed by the spherical harmonics because of the exact transformation be¬ 
tween them. In practice, however, the needlet coefficients are shrunk to zero to achieve a sparse 
representation and remove the artifacts in the function reconstructed by the spherical harmonics. 

For illustration, we apply a naive percentage-based hard thresholding method to the needlet 
coefficients at the highest frequency j — 4. The 95 percent of the needlet coefficients with smaller 
magnitudes are shrunk to zero. Figure 6 shows the function reconstructed by the spherical needlets 
after thresholding. We can see that there is no such spurious global oscillation in the reconstructed 
function because the corresponding needlet coefficients have been shrunk to zero. The values of 


these needlet coefficients are essentially artifacts that induced by the spherical harmonics. Moreover, 
thanks to the spatial localization of the spherical needlets, the thresholding of the needlet coefficients 
does not affect the reconstruction of the function within its support. More sophisticated methods 
are discussed in Scott (2011). 



Figure 5: Function reconstructed by the spherical harmonics with /max = 16. 



Figure 6: Function reconstructed by the spherical needlets after hard thresholding. 
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Figure 7: Estimated probability density functions of the needlet coefficients at frequencies j — 
2,3,4. 

Appendix: Fast Inverse Spherical Harmonic Transform 

We are interested in the following inverse spherical harmonic transform 

len I 

E E m iijk), 

l=lst m=—l 

where is one of the points on the HEALPix grid. Suppose is the p-th point (0 < p < — 1) 

on the r-th ring (1 < r < Wing)- In spherical coordinates, = {Or, (l>rp)^ where (prp = 0rO+27rp/n^. 
Using the fact that ai^m — {~^)^^lm nnd imi we have 

len I 

= I + II + II, 

l=lst m=—l 


where 


and 


I —1st 


len I 


II — ^ ^ ^ ^ ^lrn^lm{^jk)' 

l=lQt m=l 


Define Pim as the normalized associated Legendre polynomial with subscripts I and m, 

- . . I 2 I + 1 {I — m)\ ^ . 

Plm{x) = \ —- 77 -- {-yPlm{x)> 


47r (/ + m)\ 
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Then 


and 


I = ^ aioPio{cos9r), 

l = lst 


11 = 


X] X] (llmPlni{cOSer)exp{im(f)rp) 


l=lst m=l 


Interchanging the order of summation in II, we have 




len 

/=max{m,/st} 


len 

= X] fer exp(im^rp). 
m=l 


exp{im(l)rp) 


Let m = UrS + t. Then 


rir-l ( ^ 

n = ^ ^ qnrS+t,r exp [i{nrS + t)(l)ro\ > w~P 

t=0 I s J 


rir — l 


= E 


TtrW 


-pt 


t=0 


where w = exp{—27vi/rir). Thus, II can be obtained by the inverse fast Fourier transform. 
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